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Abs tract 
This  paper  presents  a  class  of  new  explicit  and  implicit  second 
order  accurate  finite  difference  schemes  for  the  computation  of  weak 
solutions  of  hyperbolic  conservation  laws.  These  highly  nonlinear 
schemes  are  obtained  by  applying  a  nonoscillatory  first  order  accurate 
scheme  to  an  appropriately  modified  flux.  The  so  derived  second  order 
accurate  schemes  achieve  high  resolution,  while  retaining  the 
robustness  of   the   original   first   order  accurate    scheme. 

"1980     Supplemented     Mathematical      Subject    Classification   65P05,    35L65, 
7  6L0  5. " 
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1.    Introduct  Ion 

In    this   paper,   we  consider     numerical      solutions      of      the      initial 
value    problem   for  hyperbolic   conservation    laws: 


Uj.   +  f(u)x  =  0    ,      u(x,0)    =  uqCx),    X  c   D. 


(1.1a) 


Here     u(x,t)      is     an     m-vector     of  unknowns  and    f(u),    the   flux,    is   a   C 
function   of    m  components.      Uq(x) ,    the  initial  data,    is   assumed    to  be      a 
function      of      bounded      total      variation      in      D.    We   assume   that    Uq(x)    is 
either   of    compact   support    in     D      =      (-",«>)      or      periodic      In     a      finite 
interval  D.    (1.1a)    can  also   be  written   in  a  matrix   form: 

u^   +  A(u)u^   =  0    ,    A(u)    =  f^    .  (1.1b) 

Equation      (1.1)    is  called   hyperbolic    if    all   eigenvalues    {a'^(u)}    of  A(u) 
are   real   and   the   set   of    its   right   eigenvectors    {r^(u)}    is    complete. 

It  is  well  known  that  solutions  to  (1.1)  may  develop  shocks  and 
contact  discontinuities  even  when  the  initial  data  are  smooth.  To 
allow  for  discontinuous  solutions,  one  admits  weak  solutions  which 
satisfy   the    system   (1.1)    in  the    sense   of    distribution   theory,    i.e., 

00 

J    J    [(tJ^u  +  <t)xf(u)]    dx   dt  +  J    (t)(x,0)   uo(x)    dx  =    0   ,  (1.2) 

0 

where  (})(x,t)    is    any   Cf°    test    function. 

Here,  and  elsewhere  in  this  paper,  we  leave  limits  of 
x-integration  and  summation  unspecified,  to  be  interpreted  according  to 
the  nature  of  the  initial  data  nQ(x)  and  its  domain  of  definition  D. 
Similarly,  the  test  functions  (t)(x,t)  are  also  assumed  to  have  the  same 
x-behavior  as   the    initial  data. 
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We     assume      that      the      system      (1.1)   possesses   an  entropy  function 
U(u),    defined   as    follows: 
(1)        U  satisfies 


Uufu=Fu  (l-3a) 


where   F   is  some   other    function   called   entropy   flux, 
(ii)      U   is    a  convex    function  of  u. 

Furthermore,  we  assume  that  weak  solutions  u(x,t)  of  (1.1)  that 
satisfy  the   entropy    inequality 

U(u)^   +  F(u)x  <   0  (1.3b) 

in  the   sense   that 

00 

-J    J    [(j)^U(u)-hJ,^F(u)]dx  dt  -  J    <t>(x,0)U(uo(x))dx  ^  0  (1.3c) 

0 

for      all     nonnegative      smooth      test      functions,      (})(x,t),      are      uniquely 
determined  by    their   initial  data  (see   [7]   and    [8]   for   more   details). 

Next  we  consider  finite-difference  (FD)  approximations  v(x,t)  to 
weak  solutions  of  (1.1).  These  are  piecewise-constant  functions 
defined  on  a   rectangular  grid: 

v(x,t)    =   v^    ,      (x,t)    e    ((j-  i.)h,(j+  i)h)    X    [(nT,(n+l)T);  (1.4a) 

h    is    the   spatial  mesh   size   and  t    is   the    time-step. 

The  values  v^  for  n  >^  0  are  generated  from  the  average  values  of 
the   initial  data 


1 


v5=^  J         u,(x)dx 

(  J-    y)h 


(1.4b) 


by 


vH+l    =    (s.v^).  (1.4c) 


where   S    is   some    FD   scheme. 

In    this    paper  we   consider  numerical    schemes   of    the    form 

L  .    v^"^l    =  R   .    v"  (1.4d) 

where  L  and  R  a  re  centered  finite-difference  (FD)  operators  that  depend 
on  {2SL  +  1)    and    (2r+l)    points,    respectively. 

If  Jl  =  0,  and  L  =  I,  the  identity  operator,  then  S  =  R  in  (1. 4c) 
and  the  scheme  is  explicit ;  in  this  case,  initial  data  with  compact 
support,  as  well  as  periodic  ones,  are  appropriate.  If  £  >  0,  then  the 
scheme  is  implicit,  and  S  =  L~^R.  The  definition  of  L"  ,  the  inverse 
of  L,  necessarily  involves  boundary  conditions;  to  avoid  the 
complication  of  treating  boundaries,  we  consider  only  periodic  initial 
data   for   implicit   schemes. 

Let  us  denote  the  dependence  of  the  FD  approximation  (1.4)  on  the 
grid,  assuming  t  =  0(h),  by  v^(x,t)  and  consider  the  limiting  process 
h  -»■   0   in   the   strip  0    <   t    <  T.    If    there   exists    u(x,t)    such   that 


lim  V|^(x,t)    =  u(x,t)    ,      Oj<t<T,      xeD  (1.5) 

h-»-0 


we   say   that   the    scheme    (1.4)    is  conve  rgent  ♦ 

When  designing   a  numerical   approximation,   we     have      to      make      sure 
that      the      resulting      scheme      is      convergent      and    that    its    limit   is    the 
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(unlque)  solution  of  the  problem.  We  would  also  like  to  have  an 
estimate  of  the  error  in  the  numerical  solution.  Unfortunately,  it  is 
difficult  to  directly  relate  the  paricular  form  of  the  discretization 
used  in  the  scheme  to  the  limiting  process  of  convergence.  Therefore, 
Lax's  equivalence  theorem  for  the  linear  hyperbolic  problem  (A  =  A(x,t) 
in  (1.1b))  played  a  major  role  in  the  development  of  numerical  schemes 
(see  [12]).  This  theorem  states  that  if  consistent,  then  L2-stability 
of  the  scheme  implies  its  convergence  in  the  L2-norm  to  the  solution  of 
the  partial  differential  equation.  Hence,  instead  of  checking 
convergence,  we  may  equivalent ly  check  the  consistency  and  the 
L2-stabllity  of  the  scheme.  With  the  aid  of  von  Neumann's  stability 
analysis,  the  questions  of  stability  and  consistency  can  be  readily 
expressed  as  design  principles.  Lax's  equivalence  theorem  also  shows 
that  the  order  of  accuracy  of  the  scheme  is  exactly  one  less  than  that 
of   its   truncation  error. 

It  is  well-known  that  L2-stability  is  necessary,  but  not 
sufficient ,  for  convergence  of  nonlinear  schemes.  Consequently, 
schemes  for  the  computation  of  weak  solutions  of  the  nonlinear  problem 
(1.1)  have  been  developed  according  to  the  following  guidelines:  (i) 
Design  a  scheme  so  that  its  linearized  version  is  L2-stable.  (ii)  Add 
numerical  dissipation  to  damp  oscillations  and  to  control  "nonlinear 
instabilities".  Since  the  addition  of  numerical  dissipation  brings 
about  loss  of  information,  the  designer  of  such  numerical  schemes  finds 
himself  in  a  position  where  he  has  to  compromise  accuracy  to  achieve 
stability,    or  vice   versa. 

The  goal   of   this   paper   is   to   show  a  way   to  avoid    this    "uncertainty 
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principle" ;    to  do    so   we  have   to  depart    from   the    linear    theory     and      "go 
nonlinear" . 

The  first  step  is  to  replace  the  requirement  of  linearized 
I'T-stabili  ty  by  that  of  stability  in  the  sense  of  uniformly  bounded 
total-variation  in  x,  and  to  add  the  requirement  of  consistency  with 
the  conservation  law  (1.1)  and  its  entropy  inequality  (1.3);  this 
implies  the  convergence  of  the  scheme  to  the  unique  weak  solution  of 
(1.1)    (see   [7]   and  Section   2). 

Next  we  use  this  convergence  theorem  to  design  schemes  that  are 
both  robust  and  accurate.  The  construction  guidelines  are  as  follows: 
(i)  Design  a  scheme  that  is  consistent  with  the  conservation  law  and 
its  entropy  inequality  and  is  total-variation  stable.  (ii)  Extract 
from  it  as  much  numerical  dissipation  as  possible  subject  to  the 
requirements    in   (i). 

The  resulting  schemes  are  second  order  accurate  (in  the  sense  of 
truncation  error  analysis)  and  genuinely  nonlinear,  i.e.,  they  are 
nonlinear  even  in   the    constant    coefficient    case. 

Numerical  experiments  reported  in  [5],  [15]-[17]  indicate  that 
this  class  of  schemes  is  superior  to  that  of  artificial  viscosity 
methods. 

The  organization  of   the   rest    of    this   paper   is   as   follows: 

(2)  Convergence   theorem. 

(3)  Basic   theory  of   total-variation-diminishing    (TVD)    schemes. 

(4)  Second   order  accurate   TVD  schemes. 

(5)  Extension   to   systems   of    conservation    laws. 

(6)  Application   to   steady  state   calculations. 
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2.    Convergence    theorem 

In      this      section     we      describe    the   convergence    theorem  of    [7]    and 
outline    its  proof. 

First   we   define    the   notions    of    consistency      and      stability      to     be 
used    in   this    theorem. 

We     say      tliat      a      finite-difference      scheme    is    consistent  with  the 
conservation  law   (1.1)    if  it    can   be    written    in   the    form 

^n+1   =   ^n   _  X(f  .+i/2-fj-l/2).  ^    =   ^ /^  (2.1a) 

where 
Vl/2    =   f(v5_,+i....,v5+,;v5ll+i,....v^:j).  (2.1b) 

the  numerical  flux,  is  a  Lipschitz  continuous  function  that  satisfies 

f (u ,  . . .  ,u;u, . .  .,u)  =  f(u)  ;  (2.1c) 

f(u)    is    the    flux  of    the   system  of    conservation    laws    (1.1). 

Similarly,  we  say  that  the  numerical  scheme  S  (1.4c)  is  consistent 
with  the  entropy  inequality  (1.3)  of  the  system  of  conservation  laws 
(1.1)  if  there  exists  a  numerical  entropy  flux  f-;+W2  that  is 
consistent    with  F(u)    (in  the    same  sense  as    (2.1c)),    such   that 

U(v^+1)  ^  U(v5)    -  A(Fj+l/2-  Fj-1/2)      •  ^^.2) 

We      turn     now     to      consider  a   sequence   of   numerical   approximations 
V,  (x,t)  ,    xe    D,    0_<t   _<T,    with  h  ->■   0  and  x    =  0(h) .      We      say      that      the 
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numerlcal  scheme  is  total-variation  (TV)  stable  if  the  total  variation 
in  X  of  v^(x,t), 

TV(v^(x,t))    =1     Iv^+i-v^l    .  (2. A) 

J 

is  uniformly   bounded    in   t   and   h;   here  n    is    the    integer   part    of    t /t . 

Theorem  2.  1  (Convergence):  Let  us  assume  that  the  scheme  (1.4)  is 
consistent  with  the  conservation  law  (2.1)  and  its  entropy  iequality 
(2.2).  If  the  resulting  numerical  approximation  is  TV  stable,  then  the 
scheme  is  convergent  and  its  limit  is  the  unique  weak  solution  of  (1.1) 
that   satisfies   the  entropy   inequality    (1.3). 

In  [7]  this  theorem  is  stated  and  proved  for  explicit  schemes 
defined  on  an  irregular  moving  mesh.  The  initial  data  considered  there 
are  of  bounded  functions  in  (-«',<=°)  that  have  bounded  total  variation. 
The  limit  (1.5)  in  this  theorem  is  in  the  sense  of  bounded,  Lj'-"- 
conve  rgence . 

In  the  present  paper  we  consider  both  explicit  and  implicit 
schemes.  To        avoid      treatment      of      boundary      conditions      we      assume 

periodicity,  unless  otherwise  stated.  Examination  of  the  proof  of  the 
convergence  theorem  in  [7]  shews  that  it  carries  over  to  the  present 
situation  with    only  minor   changes. 

LemiiH  2.2.  If  a  FD  scheme  is  consistent  with  the  conservation  law 
(2.1)  and  is  total -variation  stable  (2.3),  then  it  is  also  stable  in 
the   norm 


lv|    =    TV(v)   +   maxj    lv-j|  (2.5) 
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Proof :  We  prove  this  lemma  for  periodic  boundary  conditions. 
Since  the  assumptions  of  this  lemma  hold  for  each  of  the  components,  we 
prove  this  lemma  componentwise,  using  scalar  notation.  The  assumption 
of  total -variation  stability  means  that  there  exists  a  constant  B,  not 
depending  on   h   or  n,    such   that 

TV(v")    <  B   •    TV(v*^) 
for   all  n.    By    the    definition  of   total  variation,   we  have 


TV(v")    =  I     |v5^.i-v^|     >  max   v'j-min   vO    ; 


thus 


max   v"]   -  min     v^   <  B    •    TVCv^).  (2.6a) 

J  J 

The  combination   of   periodic  boundary   conditions      and      consistency     with 
the   conservation    law    (2.1)    implies    that    for  all   n  >^  0 


Since 


I   v5   =  I   vO    .  (2.6b) 

J  J 


rain  v"}  <  1  y    v'?   <  max   v'l    ,  (2.6c) 


where     J      is     the      number      of     mesh      points      in      a      period,    we  get    from 
(2.6a)-(2.6c)    that 
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min  v')    >   -B.  TVCv*^)    +  max  v^   >    -B.TV(v°)  +  1  Y    v^ 

J   —                                            3   —  J.J 

j                                   :  J 


ma 
J 


X   v^   <   B.TV(v^)    +  rain  v^   <      B.TV(v'^)    +  1  )]    v^; 


Hence 


max 


■    |v]|   ^B.TV(vO)    +    |i):   v^l    <B   .    TV(vO)    +  max^    \v^.\ 


and    therefore 

llv"ll    <  Bllv'^ll    ,         B  =  max(2B,l),  ' 

which    is    the   assertion  of    this    lemma. 

We  note  that  the  proof  of  Lemma  2.2  combines  the  uniform  bound  of 
the  maximal  oscillation  of  the  numerical  solution  (due  to  TV  stability) 
with  that  of  its  mean  value  (due  to  the  conservation  form)  to  conclude 
uniform  boundedness  of  the  approximation.  This  argument  remains  valid 
for  explicit  schemes  with  initial  data  of  compact  support;  it  can  also 
be   extended    to  other   boundary   conditions. 

Lemna    2.3        (Glimm) .  If        a        numerical      scheme      satisfies      the 

assumptions   of   Lemma    2.2,    then   for  all   n  and  m 

l^h(*.t„)    -^h(-.tn)lLi    <0(|t^-   tj)    . 

Proof :  See  Lemma  3.6  in  [7];  note  that  the  numerical  flux  (2.1)  is 
assumed    to  be   Lipschitz    continuous. 

We  turn  now  to  outline  the  proof  of  Theorem  2.1.  The  first  part 
of    it    follows   the    same   reasoning  as  Glimm' s   convergence    proof   in    [1]. 

Proof  of  the  convergence  theorem:  By  Lemma  2.2,  the  FD  scheme  is 
uniformly      bounded      in      the     maximum  norm.      Since    the    functions    v,  (x,t) 
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have  uniformly  bounded  total  variation  in  x  and  are  unlfoirmly  bounded, 
by  compactness  we  deduce  that  for  any  fixed  t,  a  subsequence  converges 
in  the  L^-norra.  By  a  diagonal  argument,  we  can  construct  a  subsequence 
that  is  Lj-convergent  for,  say,  all  rational  values  of  t.  Using  Lemma 
2,3  we    conclude  that    this   subsequence   converges   in   the   norm 


max   II  v(«  .t)!!^ 
t  ^ 


to  a  limit  that  we  denote  by  u(x,t).  Certainly,  u(x,t)  is  uniformly 
bounded  and  has  a  uniformly  bounded  total  variation  in  x.  By  Theorem 
1.1  of  [8]  (which  is  a  straightforward  extension  of  the  Lax-Wendroff 
theorem  in  [11])  u(x,t)  is  a  weak,  solution  of  (1.1)  that  satisfies  its 
entropy  inequality    (1.3). 

Finally,  we  claim  that  not  only  a  subsequence  of  v^(x,t)  but  the 
whole  sequence  converges  to  u(x,t).  For  if  not,  then  by  cmpactness  we 
could  select  a  subsequence  that  tends  to  a  limit  u  ?^  u.  Now  u  and  u  are 
both  weak  solutions  of  (1.1)  that  satisfy  its  entropy  inequality  (1.3); 
furthermore,  they  both  have  the  same  prescribed  initial  data  Uq(x). 
Therefore,  by  the  uniqueness  of  the  initial  value  problem  assumed  in 
Section  1,  we  conclude  that  u  e  u.  This  proves  that  the  scheme  is 
convergent . 

Note  that,  because  of  the  indirect  nature  of  the  convergence 
proof,    we  have   no   error   estimates. 
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3.    Basic    theory  of    TVD   schemes. 

In  this  section  we  consider  the  scalar  conservation  law  (m  =  1  in 
(l.D) 

u^   +  f(u)^  E    ut   +  a(u)ux  =    0  .  (3.1) 

Equation  (3.1)  states  that  the  solution  is  constant  along  the 
characteristic  lines  dx/dt  =  a(u).  Hence,  if  the  initial  data  are 
smooth,  the  total  variation  in  x  of  the  solution  remains  constant  in 
time,  as  long  as  the  characteristic  lines  do  not  intersect.  Even  when 
a  shock  forms  due  to  convergence  of  characteristics,  the  solution  at 
each  point  (x,t)  remains  well-defined  in  terms  of  a  backward  drawn 
characteristic  line,  connecting  it  to  the  initial  data.  Since  the 
initial  data  in  the  domain  of  dependence  of  the  shock  curve  are 
dissipated  into  it,  the  formation  of  a  shock  may  actually  decrease  the 
total   variation   in   x    (see  Fig.       1). 

In  view  of  this  behavior  of  the  total-variation  we  are  interested 
in  designing  numerical  approximations  to  the  scalar  conservation  law 
(3.1)    that   are    total -variation   diminishing    (TVD). 

In  this  section  we  review  the  theory  that  has  been  developed  in 
[5]    for   explicit  TVdI    schemes,    and   extend    it   to   implicit    ones. 

^ In  [5]  we  use  the  term  TVNI  (total  variation  non-increasing)  instead 
of    TVD. 
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Flgure   1 


We    define    the   total   variation  of  a    mesh  function  v  to  be 


TV(v)   =  I       |A^i/2v| 


(3.2) 


where  A  j,.w2V  =  V:j^.^-v^.      (We   use   the    general  notation   convention 


^j+l/Z^  =  ^j+1    "  ^j 


(3.3) 


for  any  mesh  function  b.)  We  say  that  a  numerical  scheme 


,n+l  =  s  •  v" 


(3.4a) 


is  TVD  if 


TV(v"+l)  <  TV(v")  ; 


(3.4b) 
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here  v"   denotes    an  approximation   to   the    solution  of    (3.1)   at      t      =     ni , 
where  t    is    the    time-step. 

In      this      paper  we    consider   the   one-parameter    family   of   schemes   in 
conservation    form 

L  .    v""^^    =  R   .    v'^  (3.5a) 

where   L   and  R  are   the    following    finite-difference   operators: 

(L.v)  .   =  Vj  +  nX(?j+i/2   -  ^j-1/2)  (3.5b) 

(R.v)j   =  Vj    -   (l^)A(?j+i/2    -^j-1/2)    •  (3.5c) 

Here  f  is  a  numerical  flux  that  is  consistent  with  f(u)  in  (3.1).  The 
scheme  (3.5)  can  be  written  in  the  form  (2.1)  with  f  =  (l~n)f^j.i/2  + 
^^i+1/2'  ^^^  therefore  is  consistent  with  the  conservation  law  in  the 
sense  of  (2.1).  The  scheme  (3.5)  with  n  =  0  is  the  explicit 
forward-Euler  scheme;  for  n  >  0  (3.5)  is  implicit,  thus  S  =  L~^R  in 
(3.4a).  The  scheme  (3.5)  with  n  =  1  is  the  implicit  backward-Euler 
scheme. 

LemnH    3.1.      If  R    is  a   total -variation   diminishing    (TVD)      operator, 
i.e.  , 

TV(R  .    v)  _<  TV(v)  (3.6a) 

and  L   is  a   total-variation   increasing    (TVI)   operator,    i.e., 

TV(L   •    v)    >  TV(v),  (3.6b) 
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then  the  numerical   scheme    (3.Aa),    with  S  =  L~^R,    Is  TVD. 

Proof  ;    TVCv^+b  _<  TV(L.  v"+l)    =  TV(R.v")  _<  TV(v"). 

In    this   section  we  consider  3-polnt   schemes.    I.e.,    f 4+1/2   ^^  (3-5) 
is  of   the   form 

f"j+l/2   =  f(Vj,Vj+i)  (3.7a) 

where 

7(u,u)    =  f(u)    .  (3.7b) 

Assuming    the   numerical    flux  in   (3.7)    to  be      Lipschitz      continuous, 
we   can  write  It    in   the   general   form 


fj+1/2  =   f(vj,Vj+i)    =j  (fj+fj+l)    -  2X  ''j+l/2^J+l/2''  ^^-^^ 


where     q4+W2      =     q(v^,v^+p      is     some  bounded    function  and  ^i^\/2^   ^^ 

(3.3).      The   notation   f^   =   f(v^)   is  used  throughout   this   paper. 

Another   general    form   to  describe  a  Lipschitz    continuous      numerical 
flux    in  (3.7)    is 


fj+1/2   =    f(^j.Vj+i)    =  fj   -^ct^.i/2Aj+i/2V 


"  ^j+1  -  x   cT^l/2^  j4-l/2^»  (3-^^) 


where  C.^W2    =  C    (v.jV^+j^)   are   some   bounded  functions.      Comparing    (3.8) 
with  (3.9a),    we   get 
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^j+1/2    =   2    ^''j+l/Z  ^   ^  j+1/2)  (3.9b) 


where  "^^i/o   i^    the  mean-value   local   CFL   number 


^j+1/2    =  ^   3j+i/2  (3.10a) 

^j+1/2   =  ^  j+l/2f    /   ^jfl/2^   •  (3.10b) 


Using  (3.9a)  for  the  numerical  fluxes  in  (3.5a)  and  (3.5b),  we  see 
that  the  finite-difference  operators  L  and  R  can  be  expressed  in  the 
form 

(L.v)j   =  Vj    -  n(ct^i/2Aj+i/2V  -  Cj_i/2Aj_i/2v),  (3.11a) 

(R.v)j   =  Vj   +  (l-n)(Cj+i/2Aj+i/2V-C]_i/2Aj_i/2v).  (3.11b) 

Next  we  formulate  sufficient  conditions  for  a  finite  difference 
operator  to  be   TVD   (3.6a)    or  TVI    (3.6b). 

LemnH    3.2.      Let   Z   be  a    finite    difference   operator   of   the    form 

(Z.v)j   =  Vj   +  ct^i/2  Aj+i/2V  -  CT_^/2  Aj_i/2V  (3.12a) 

where  A^^]^/2V    is  (3.3). 

(a)  If    for  all    j 

Cj+1/22  0.         ct^l/2  +  CT+1/2   <    1  (3.12b) 

then  Z   is  TVD. 

(b)  If    for  all    j 
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1  ^j+1/2  1 


(3.12c) 


then  Z   is  TVI. 

Proof :      Denote   w.   =   (Z«v)::   and   subtract   (3.6a)   at    j   from   (3.6a)   at 
j+1    to  get 

Aj+i/2W  =    (l-Cj+i/2-C]+i/2Mj+l/2^  +  ct^3/2Aj+3/2V  +  C~yi/2^ j-U^'''    (3.13a) 

The      Inequalities     (3.12b)      are  exactly        the        conditions        for        the 

coefficients      of     ^t+]^+i/2'^*      ^  ~      ^»~^      •'■"   (3.13a)    to   be   nonnegative. 

Hence    taking    the   absolute    value  of      (3.13a)      and     using      the      triangle 
inequality  we  get 

l^j+l/2w|   i  (l-Cj+l/2-Cj+l/2)l^j+l/2v|    +  ct^3/2l^j+3/2l    "^^^  j-1/2  1^  j-l/2v|    ' 

(3.13b) 

summing      (3.13b)      over      j,      assuming      v      to     be      periodic  or   of   compact 
support,    we   obtain 


TV(Z.v)    =1    |Aj+i/2w|    <l    (l^j+i/2-C]+l/2)|Aj+l/2v| 
j  j 

+  1    Ct^3/2l^j+l/2v|    +1    CT_^/2lAj-l/2^l 
j  j 


this   completes   the   proof    of   part   (a). 

To    prove   part  (b)  we    rewrite  (3.13a)   as 


J 


-19- 
Aj+i/2W+    (-Cj+3/2)A  j+3/2^  +    ("^  j-1/2)^  j-1/2^  =    <^l-Cj+l/2   -  ^]-l/2^^  j^-l/2''    (3.13c) 

and  observe  that  (3.6c)  is  the  condition  for  the  coefficients  of 
^  T+ir+]/2^  ^  "^  0,±1  in  (3.13c)  to  be  nonnegative.  Hence  taking  the 
absolute    value  of    (3.13c)    and  using    the    triangle   inequality,   we  get 


l^j+l/2w|    -  ct^3/2|A^3/2v|    -C    j-1/2  |A  j-i/2v| 

2(1   -   ct+1/2   -  C]+i/2)|A^l/2v|    .  (3.13d) 


Summing    (3.13d)   over    j,    we  obtain 

TV(Z.v)    -I    (ct+i/2+C]+i/2)|Aj+i/2v| 
2 

=  1    |Aj+i/2w|    -I    Cj+3/2lA-^3/2v|    "I    ^'j-l/l    l^j-1/2^1 
J  J  J 

ll    (l^j+l/2-Cj+l/2)    l^j+1/2^1 
j 

=  TV(v)    -I    (ct+1/2   +CJ+1/2)     |Aj+i/2  v|    . 

j 

Comparing   the   two    extreme   parts   of    the  above   relation  we   conclude   that 

TV(Z«v)  2  TV(v)    ; 

this   completes    the    proof   of    Lemma    3.2. 

Next  we  use   Lemma   3.2   to   formulate   a    sufficient    condition    for      the 
schemes   (3.5)    to   be   TVD. 

Theorem   3.3.      If    qj+1/2    i"   (3.8)    satisfies 


^  (3.14) 


'J+I/2I   iqj+1/2  <T^  . 
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where     '^4+1/2   i^    the  rnean   value    CFL-number   (3.10)    then    the    scheme    (3.5) 
with  (3.8)    is   TVD. 

Proof:    Comparing    (3.11)    with  (3.12a)   we    see    that    both      the      finite 
difference   operators   L  and  R  can  be    written    in    the    form   (3.12a)  with 


C-    =    (l-n)C^      for  R  (3.15a) 

C-   =    -  n    C-        for  L  (3.15b) 


Using  Lemma  3.2  we  conclude  from  (3.9b),  (3.14)  and  (3.15)  that 
the  operator  L  in  (3.5)  is  TVI,  while  the  operator  R  is  TVD. 
Therefore,    it   follows   from  Lemma    3.1    that  the    scheme    (3.5)    is   TVD. 

In  the  following  we  consider  the  class  of  schemes  (3.5)  and  (3.8) 
with 

qj+1/2   =   Q(^j+l/2;^j+l/2)  (3.16a) 

where   Q(x;e)    is    either 

[  e         ,      |x|    <  e 
Q(x;e)    =       J  (3.16b) 

[  |x|     ,       |x|    >   e 
or  9 

f  i-  (—  +  e)    ,       |x|    <  e 

Q(x;£)    =  2      e  .  (3.16c) 

[  |x|  ,       |x|    >   e 

here  e  i^.W2    -  ^(^i'^i+l)  2.  ^*    ^^^   that 

Q(x;0)    =    |x|      and      Q(x;e)    >e    .  (3.16d) 
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Whenever  we    are    not    concerned   with   the   c    dependence   we    shall     use     Q     = 

Q(x).         For     e    >    0  both  (3.16b)   and    (3.16c)   are   positive   approximations 

to    |x|;    Q(x)    in   (3.16c)    is  continuously   dlf ferentlable,    but   it    is      only 

Llpschitz      continuous      in      (3.16b).        If    0  <   e    <  in   (3. 16b)-(3. 16c) 

~  1-n 

then     1^+1/2      i^i      (3.16a)      satisfies      the      inequality      (3.14)      for     all 

1"^  4+1/9  I     <  -, .         Therefore,      it      follows   immediately  from  Theorem   3.3 

jT '■I  ^     —   1-ri 

that 

Corollary  3.4.  The  scheme  (3.5)  and  (3.8)  with  q-:^.W2  defined  by 
(3.16)    is    TVD  under    the    CFL-like   restriction 

X    max    |a^i/2l    <  JZ^  •  (^-l^) 

Observe  that  the  backward -Euler  implicit  scheme,  n  =  1  in  (3.5), 
is  unconditionally  TVD,  while  the  f orward-Euler  explicit  scheme,  r,  =  0, 
is  TVD  under   the  CFL  condition  of   1. 

In  [5]  we  show  that  the  set  of  TVD  schemes  contains  the  set  of 
monotone  schemes.  However,  unlike  monotone  schemes  (see  [6]),  the  fact 
that  a  scheme  is  TVD  does  not  automatically  imply  that  it  is  consistent 
with  the  entropy  inequality  (2.2).  For  example,  the  first  order 
accurate  upstream  differencing  (3.5)  and  (3.8)  with  q^+1/2  =  1^1+1/21 
(i.e.  e  =  0  in  (3. 1  6b)-(3. 1  6c)  )  admits  stationary  "expansion  shocks" 
(see   [8]  and    [9]). 
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In  [3]  and  [9,  Appendix  A]  we  represent  the  numerical  flux  (3.8) 
with  (3.16)  as  the  flux  across  x  =  0  In  an  approximate  solution  to  the 
Riemann  IVP 


UqCx)    = 


^j 


X  <   0 


Vj+i    ,      X  >  0 


(3.18a) 


The       approximate        solution      to      the      Riemann     problem     w(x/t  ;v .  ,v^j.i ) 
corresponding    to  e    =    0,    i.e.      Q(x)    =    |x|,    is 


w(x/t;Vj,Vj+i)    = 


V.        ,      x/t    <  a^i/2 


(3.18b) 


^j+1  '  '^/t  >  aj+1/2 
where  a^+W2  is  (3.10b).  The  fact  that  (3.18b)  approximates  the 
solution  to  (3.18a)  by  a  shock,  regardless  of  the  entropy  condition,  is 
responsible  for  the  stationary  "ejqjansion  shocks"  in  the  upstream 
differencing  schemes  mentioned  above.  Therefore,  we  consider  the 
following  modification  of    (3.18b), 


(Vj  ,      x/t    <  aj+i/2-6 


w(x/t;vj,vj+i)=' 


Vj+l/2(^/t)'      aj4-l/2-'5    <  x/t    <   a^Y/2-^ 
yj+1  .      aj+l/2+'5    <  x/t 


(3.18c) 


(3.16b)   corresponds    to  taking  v.^2./2(x/t)   to  be   the   constant   state 
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Vj+i/2(x/t)    =   ^  (v^  +  vj+i)    ,  (3.18cl) 


while    (3.16c)    corresponds   to    the    linear    function 


^j+l/2(x/t)    =  J  (Vj+  Vj+i)    +  2^ (Vj+i-  vj)    ;  (3.18e) 


6   has    the   dimension   of   velocity  and 

e    =  X   6  (3.18f) 

in   (3.16). 

In  [9,  Appendix  A]  we  show  that  if  6  in  (3.18c)  is  chosen  so  that 
the  fan  |x/t  -  Si+i/ol  '^  ^  covers  the  domain  of  influence  of  the 
initial  discontinuity  (3.18a),  then  the  approximate  solution  in 
(3. 1 8c)-(3. 1 8d)  is  consistent  with  the  entropy  inequality  of  the 
Riemann   problem. 

Using  tlxe  representation  of  the  forward-Euler  scheme,  n  =  0  in 
(3.5),  (3.8)  and  (3.16b),  as  a  Godunov-type  scheme  corresponding  to  the 
Riemann  solver  (3.18c)  (see  [8]  and  [9])  we  deduce  that  it  is 
consistent  with  the  entropy  inequality  (2.2).  Since  this  scheme  is 
also  cnsistent  with  the  conservation  law  (2.1)  and  it  is  TVD  under  the 
CFL-restriction   (3.12)   of    1,    we   conclude  by  Theorem  2.1    that 

Corollary  3.  5.  The  explicit  forward-Euler  scheme,  n  =  0  in  (3.5), 
(3.8)  and  (3. 1 6a)-(3. 1 6b) ,  is  convergent  and  its  limit  is  the  unique 
weak   solution  of    (3.1). 
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We     note      that      the   modification  of    |x|    in   (3.16b)is  equivalent    to 

adding  an  entropy   viscosity   term  with  a   coefficient   e    to  the     upstream 

differencing      scheme.        Choosing     e      =     ^  ^-(-1/2   locally    by  the   domain  of 
influence   condition  mentioned  above,   we  get: 

(i)  ej+1/2    =0(l^j+l   -  ^jl)  (3.19a) 

(ii)  ^rl-1/2    =   0     if      (3.18a)   is  a  shock.  (3.19b) 

Hence  the  added  entropy  viscosity  term  is  0(|Av|  )  and,  unlike  the 
standard  artificial  viscosity  terras,  it  vanishes  in  the  neighborhood  of 
shocks. 

The  schemes  in  the  class  (3.5)  with  (3.8)  are  all  first  order 
accurate  in  both  time  and  space,  except  for  n  =  — ,  where  the  scheme  is 
second  order  accurate    in   time. 

In    the   next   section  we      present      a      rather      general      technique      to 

convert      a     3-point      first    order  accurate   TVD  scheme   of    the    form    (3.5), 

(3.8)    into  a  5-point    second   order  accurate    (in  both   time   and   space,      or 

just      space)   TVD  scheme   of   the   same   generic    form.      This  second   stage   of 

our   design  of  high-resolution  TVD  schemes    capitalizes   on   the    fact      that 

the   exact  solution   to    (3.1)    is   TVD  due   to    the    phenomenon   of   propagation 

along   characteristics,    and    is   independent    of    the   particular   form  of    the 

flux      f(u)      in      (3.1).        Similarly,      the      original    first    order  accurate 

scheme    is  TVD  subject    only  to   the   CFL    restriction    (3.17),    independently 

of      the     particular      form     of      the      flux.      Thus  to  achieve   second  order 

accuracy   while   retaining    the   TVD  property,      we      use      the      original     TVD 

scheme   with  an  appropriately   modified    flux   f   +  _  g.    The    requirements    on 

A 
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g  are:  (1)  The  characteristic  speed  Ag/Av  should  be  uniformly  bounded. 
(ii)  The  modified  .  scheme  should  be  second  order  accurate  almost 
everywhere.  The  details  of  the  construction  of  such  a  g  are  given  in 
the   next   section. 

We  would  like  to  point  out  that  the  resulting  second  order 
accurate  scheme  has  to  be  genuinely  nonlinear,  i.e.  nonlinear  even  in 
the  constant  coefficient  case,  since  linear  second  order  accurate  TVD 
schemes    do   not    exist    (see    [5]). 
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4,    Second  Order  Accurate  TVD  Schanes 

In  this  section  we  convert  the  first  order  accurate  TVD  schemes  of 
Corollary  3. A  into  second  order  accurate  ones  by  using  a  technique 
developed  in    [ 5] . 

The  basic  idea  is  to  use  the  first  order  accurate  scheme  with  an 
appropriately  modified    flux   f   +  _ g,    i.e. 

A 


^^j+l/2   =  2    ^^j"^  ^j+l^   +2"  ^®J^J+1^"  2  ^^''j+1/2  +  ^j+l/2>  ^j+1/2^    ^^'l^) 


where 

"^^4-1/2   =  ^j+1/28    /  ^j+1/2  ^    »  i^.lh) 

thus   preserving    the   TVD  property    of   the    scheme. 

The  construction  of  g  to  be  used  in  (4.1)  is  done  in  two  steps: 
First  we  use  truncation  error  analysis  to  find  the  Taylor  expansion  of 
g  for  which  (3.5)  with  (4.1)  is  a  second  order  accurate  approximation 
to  (3.1).  Then  we  apply  a  regularization  procedure  to  the  point  values 
of   g   to   ensure   that  Y  i+W2   in   (4.1b)    is    bounded. 

The  first   step   is   summarized    in   the    following   lemma: 
Lemna    4.1.      If  g   in   (4.1)   has  the  Taylor  expansion 

g   =  h  a(v)   v^  +  0(h2)  (4.2a) 

where  v   =  Xa    and 

a(v)    =1q(v)   +  (n    -i)  v2  (4.2b) 
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then  (3.5)    is  a   second  order  accurate   approximation    to    (3.1). 

Proof:      The      function  Q(x)(3.16)    is  at    least   Lipschitz   continuous, 
therefore 

|Q(v^l/2+^j+l/2)    -Q(^j+l/2)l     l^j+1/2^1 
<  const.     lTj+i/2l     |Aj+l/2^l 


=   const.    |g.^^   -   gJ      =     0(h    ) 


Hence  the   Taylor  expansion    to  0(h2)   around     Xj^^/2      °^      (4. la)      can     be 
written   as 

xT.^1/2    =[Xf  +  g   -^Q(v)Vx]j+l/2   +0(^2)  (4.3a) 

using    (4.2)  we  get 
xTj^l/2    =[Xf  +h(n    -|)v2vJj^i/2   +0(^')-  ^'-^^^ 

It    follows  therefore    that 

A(?5^,/2   -  f  j-1/2)    =   t^fx  +    (n-  ^)   T2(a2v^)J5  +  0(h3) 
X(F^:l/2-T5ll/2)    =X(F^^i/2-   ^5-1/2)    ^x2(f,,)5-.0(h3)    . 

where      we  have  assumed    the   0(h2)    terras    in   the  above   expansions   to  be  at 
least   Lipschitz    continuous  and   t    =  0(h). 

Using   these   expansions    in  (1.8a)   we    get   by   rearranging    terras    that 

v^+l  =  v5  -   [Tf^  +   (n-  i)  T2(a2v^)J]  -  n  x2(f^^)n  +  0(h3) 


-28- 
=    [v   -  T    f^  +  it2(a2v^)^]'l  +  0(h3) 

which   completes    the    proof   of    the    lemma. 

We   remark   that    if   a  (v )   in   (4.2b)    is    taken   to  be 


o(v)    =1  Q(v)  (4.4a) 


then  (4.3b)    becomes 


A    fj+1/2    =A    fj+i/2   +0(h'^) 


2^    . 


(4.4b) 


which  shows  that  the  modified  scheme  is  second  order  accurate  in  space; 
this   observation  will   be   used   when   dealing   with  steady  state    solutions. 

The  essential  part  of  the  second  step  is  given  in  the  following 
lemma: 

LemnH    4.2.      Let   g   be  defined  by 

gj   =    s«max{  O,rain(o  j+j/2  l'^j+1/2^1  »^*<' j-1/2  ^j-1/2'^^}  (^.53) 

where 

s  =  sgn(A  j^^/2  v)  (4.b) 

and  o^W2    =  ^^^i+1/2^    ^^   (4.2b).      Then 

(a)  g  =  h   a(v)v^  +  0(h2)  (4.6) 

(b)  ITj+1/21    =    |Aj+i/2   g    /  ^j+1/2^1   i<'(^j+l/2>  ^^-^^ 
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(Note  that   a  (v  )    >  0  for    |v  |    <  _ ), 

Proof :    (a)    See    [5],    or  verify    directly, 
(b)   Note  that    the    sequence    g.    in     (4.5a)      cannot      change      sign     without 
vanishing  at  the    transition   point.      Consequently 


|gj+l"  gjl   1^^    (ISjI.    ISj+ll  ) 

<  max    [min   (a  ^1/2  |A  3+1/2^!  ,    a  j-1/2  1^  j-1/2^1  ^'   ™i"(a  j+3/2  1^  j+3/2v| 
°  j+l/2lAj+l/2  v|)]  i  Oj+1/2    |Aj+i/2   v|    . 

Comparing  the  two  extreme  parts  of  the  above  inequality,  we  get 
(4.7).      Vfence    the   characteristic    speed  Y4+W2   is   uniformly  bounded. 

The  most  general  g  in  (4.1)  to  satisfy  the  requirements  of  our 
construction   is    of    the    form 

g  =   g  +  g  (4.8a) 

where    g   is   (4.5)    and  g    is   subject    to   the    requirements 

(1)  g   =  0(h2)  (4.8b) 

(ii)  '<'t+-1/2    ^  ^  i+l/2S''^  i+1/2'^    ^^   uniformly   bounded.  (4.8c) 

We  summarize  the  results  of  our  construction  technique  in  the 
following  theorem. 

Theorem  4.3.  The  scheme  (3.5)  with  (4.1)  where  g  is  given  by 
(4.5)  and  (4.8)  is  a  second  order  accurate  approximation  to  (3.1)  and 
is  TVD  for  all  0  <  n    <   1  under   the  CFL-like   restriction 
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raaxj    |vj^.i/2   +  Yj+i/2l   i  yir  •  (^-^a) 


Proof:  It  follows  from  (4.8b)  and  (4.6)  that  g  in  (4.8a)  satisfies 
(4.2).  In  [5]  we  show  that  the  coefficients  in  the  0(h2)  term  in  (4.2) 
may  be  discontinuous,  but  only  at  a  finite  number  of  points.  Hence  by 
Lemma  4.1    the    scheme   is   second  order  accurate. 

That  the  scheme  is  TVD  follows  from  its  being  the  same  TVD  scheme 
of   Corollary    3.4,    except   that    it    is   applied    to   the    flux   f   +  _  g. 

A 

Exactly  as  in  (3.9),  the  numerical  flux  (4.1)  can  be  rewritten  in 
the    form 

^    Vl/2    =  ^    ^  j  ^  8j    "  ^j+1/2  ^j+1/2   ^   =  ^^j+1"^  8j+l"  Cj+l/2^j+l/2^(^-^b) 
where 


Therefore  we  conclude  by  the  same  way  as  in  Corollary  3.4  that  the 
modified  scheme  is  TVD  under  the  modified  CFL  restriction  (4.9a);  this 
completes    the   proof    of  Theorem  4.3. 

Unless  otherwise  stated  we  take  in  (4.8)  g  H  0;  hence,  from  this 
point   on,    g  =    g   (4.5)   and  y   =   y. 

We    remark    that   by   (4.7) 

|v    +  yI  i    |v|    +    \y\    <    |v|    +  o(v)    ; 

hence  the  CFL  restriction  on  |v+y |  (4.9)  can  be  expressed  in  terms  of  a 
CFL  restriction   on    |v  | : 
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maxj    |v^^/2l  i  V  (4.10a) 

where      vj      is      the    largest   value   of   x  for  which   the    following    inequality 
holds 


|x|    +  0(x)    <  J_   ;  (4.10b) 

—  i-n 


a(x)   is   (4.2b). 

We  observe  that  in  the  two  important  cases  of  n  =  0  and  n  =  1,  the 
restriction  (4.10)  is  identical  to  the  original  one  (3.17),  i.e.,  p  =  1 
for  n  =  0  (see  [5])  and  y  =  ^  for  n  =  1. 

Next  we  discuss  a  particular  property  of  the  numerical  flux  (4.1) 
with  (4.5).   The  dependence  of  g.  on  v  in  (4.5)  is  of  the  foi 


)rra 


gj   =   g(vj_^,Vj,Vj+i)    .  (4.11a) 

^1+1/2      ^^   (4.1)    depends   on   g^  and   g -i-i-i    ,    therefore    its    dependence   on   v 
is   of   the    form 

^j+1/2   =  f(^j-l.Vj.Vj+i,Vj+2)    •  (4.11b) 

Since  (3.5)  depends  on  f^_i/2  arid  f-i+i/2  >  ^^^  modified  scheme 
becomes  a    5 -point   scheme. 

The  consistency  relation  of  a  numerical  flux  of  a  genuinely 
5 -point    scheme    is 

f(v,v,v,v)    =  f (v)    . 
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The  numerical  flux  of  our  second  order  accurate  TVD  scheme, 
however,  has  a  consistency  relation  which  is  similar  to  that  of  the 
original   3-point    scheme    (3.7),    in  the    following    sense: 

Property  4,4.  f-j+i/2  (4.11b),  the  numerical  flux  (4.1)  with 
(4.5),    satisfies 

T(u,v,v,w)    =  f(v)  (4.11c) 

for  all   V  arri    for  all  u  and   w. 

Proof  :  The  particular  form  (4.11c)  is  equivalent  to  setting  v.  = 
^4+^  =  V  in  (4.1)and  (4.5);  thus  Ai^.w2V  =  0.  In  this  case  g(u,v,v) 
g(v,v,w)  =  0  for  (4.5)  in  the  form  (4.11a),  or  equivalently  g.  =  gj^.^  = 
0  in  (4.1).  Since  v  ^^^ /2  (3.10)  aridYi-(-i/2  (4.1b)  are  bounded  (see 
(4.7)X,  it  follows  that  Q(v -i+W2  "'"Y-h-1/2)  ^i+1/2  ^  =  0  (independent  of 
the  particular  value  assigned  to  V4+]^/2  ^"^  'Y  i+1/2  ^°^  ^  i+1/2  ^  =  0). 
Hence  fj+1/2  =2  tf(vj)  +  f (vj+d  ]  =  f{v),  which  completes  the  proof  of 
property   4.4. 

Having  property  4.4  in  mind,  we  shall  loosely  refer  to  this 
5-point  modified  scheme  as  "essentially  a  3-point  scheme".  This 
property  is  the  reason  for  the  ease  of  handling  the  extra  numerical 
boundary  conditions,  and  for  the  tridiagonal  structure  of  the 
linearized  LHS  of   the    implicit   schemes    (to   be    discussed   later). 

Now  that  we  have  shown  that  the  modified  second  order  accurate 
scheme  is  consistent  with  the  conservation  law  (3.1),  we  conclude  by 
Theorem  4.3  and  the  first  part  of  the  proof  of  the  convergence  theorem 
that 

Corollary   4. 5.         the      second      order  accurate   TVD  scheme   of   Theorem 
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4.3  has  a    convergent    subsequence,    the    liTii  t   of   which   Is    a   vv^eak   solution 
of   (3.1). 

Numerical  experiments  reported  in  [5]  indicate  that  the  second 
order  accurate  TVD  scheme  of  this  section  (with  an  appropriate  choice 
of  e  in  (3.16),  (3.18))  is  consistent  with  the  entropy  inequality; 
however,  at  present  we  do  not  have  a  rigorous  proof  for  this  fact. 
Consequently  we  cannot  use  the  second  part  of  the  equivalence  theorem 
to   conclude   full  convergence.      We  console    ourselves   with    the    following 

Corollary  4.6.  The  second  order  accurate  TVD  scheme  of  Theorem 
4.3   is    convergent    in    the    constant    coefficient    case. 

This  follows  immediately  from  the  uniqueness  of  the  initial  value 
problem   in  the    constant    coefficient    case. 

Again,  we  would  like  to  point  out  that  the  second  order  accurate 
TVD  scheme  is  nonlinear  even  in  the  constant  coefficient  case  and 
therefore.  Corollary  4.6  cannot  be  obtained  via  a  standard  linear 
l'2~stabili ty   argument. 
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5.    Systems    of   Conservation  Laws. 

In  this  section  we  describe  how  to  extend  the  new  second  order 
accurate  scalar  schemes  of  Section  4  to  systems  of  conservation  laws 
(1.1).  Our  extension  technique  is  a  somewhat  generalized  version  of 
the   procedure   suggested  by  Roe   in    [13]. 

Let 

R(u)    =  [rku),...,r"i(u)]  (5.1a) 

be  an  raxm  matrix  the  columns  of  which  are  the  right  eigenvectors  of  the 
Jacobian  matrix  A(u)  (1.1b).  Since  the  set  of  right  eigenvectors  of 
A(u)  is  assumed  to  be  complete,  R(u)  (5.1a)  is  invertible.  We  note 
that  the  rows  £  ^u) ,  .  . .  ,£™(u)  of  R~^(u)  constitute  the  bi-orthonormal 
system  of    left-eigenvectors    of    A(u), 

£^rJ   =  6. J  (5.1b) 

and    that 

R"l  A   R  =  A    ,  ^ij   =  a^'^ij    ♦  (5.1c) 

where    {a"'"(u)}    are   the   eigenvalues   of   A(u). 

It      is      shown      in      [10]    that   the   initial  value   Rieraann  problem   for 
(1.1) 


Ut     ,        x   <    0 
Uo(x)    =        \  (5.2a) 

Ur     ,  X    >    0 

where    Ii^d-Ut  |    is    small,  has    a   unique   solution     which      consists      of      m+1 
constant   states 
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u^  =   u°,    u^,    ....  u"     =     u^  (5.2b) 

separated  by  m  centered  waves.  If  the  k-characteristlc  field  is 
genuinely  nonlinear  (a'^r'*'  ?^  0),  then  the  k-wave  separating  u*^"^  and  u'^ 
Is  either  a  shock  or  a  rarefaction,  depending  on  whether  the 
characteristic  field  is  convergent  or  divergent;  if  the  k-fleld  is 
linearly  degenerate  (a^r*^  =  0)  then  the  k-wave  is  exclusively  a  contact 
discontinuity.      The  intermediate    states    {u^}    satisfy 

uk   _  ^k-1   =   ~k    j.k(^^)    +0(|ur-ul|2)  (5.2c) 

where 

^   =  £^(ul)(ur  -  Ul)    .  (5. 2d) 

Unlike  the  scalar  case,  the  total  variation  in  x  of  the  solution 
to  the  system  of  equations  (1.1)  is  not  necessarily  a  monotonic 
decreasing  function  of  time,  and  it  may  actually  increase  at  moments  of 
interaction  between  waves. 

Ideally,  we  would  have  liked  to  extend  the  scalar  concept  of  TVD 
schemes  to  systems  through  some  diminishing  functional,  that  bounds  the 
total  variation  in  x  from  above,  i.e.  the  kind  of  functional  that  was 
constructed  by  Glimm  in  [1]  for  initial  data  of  small  total  variation. 
At  present,  we  do  not  have  such  a  monotonic  decreasing  functional,  nor 
any  other  vehicle  that  apriori  ensures  the  uniform  boundedness  of  the 
total   variation    in   x.    So,    we  do   what    we   can   do... 

In     the      following     we      describe      how  to  extend   our   new  scalar  TVD 
scheme    to   systems   of    conservation    laws    so  that    the   resulting   scheme      is 
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TVD       for        the        "locally     frozen"      constant-coefficient      system.        To 
accomplish      that      we      define      at      each      point        a        local        system        of 
characteristic      fields,      and    apply   our   scheme    "scalarly"   to  each   of    the 
fields. 

Let  v.^,  /2   denote   some   symmetric   average   of    v^   and   v^^^    ,    i.e. 

^j+1/2   =  V(Vj,Vj+i)  (5.3a) 

where 

V(u,v)    =  V(v,u)    ,  V(u,u)    =  u    ,  (5.3b) 


and   let   a^+i/2>    '^i+l/Z   ^"^  '^i+l/Z   denote   the   respective      quantities      of 
A(v.^^  ,2)'      Ijet  OL^+i/2   ^^   defined   by   the    relation 

m 

Vj+1   -   v-j   =     I      a'5+1/2   rH^i/2    ;  (5.4a) 

k=l 

it  follows   from   (5.1b)    that 

"j+1/2   =  ^j+l/2(^j+l    -^j)    •  (5.4b) 


°'i+l/2      ■'■^     exactly      the    jump  of    the   k-wave    in    the   constant-coefficient 
Rieraann  problem 

r   V .      ,    X  <  0 

Ut+A(v.+,/2)u   =    0,      u(x,0)    =    \  ;  (5.4c) 

I     ^j+1    ,      X  >  0 

it    is    an  0(|Ai^./2v|  )   approximation    to      the      size      of      the     k-wave      in 
(5.2). 

We     now     extend  the    scalar   scheme   of    Theorem  4.3   to   general    sysems 

of    conservation    laws  as   follows. 
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vf  1   +  nX(F^:j/2-  TfJy2^    =  v5   -   (1-n)   X(7^^i/2-  ?^_i/2)  (5.5a) 

where 

m 

"^j+1/2   =y  {fj+fj+1   +^     I      [gj   ^gj+l   -Q^+1/2   +Yj+l/2)«j+l/2]'^j+l/2H(5.5b) 

k=l 

g j   =  s    •    max    [0,   min(a(vH_|_^^2)  |aj+i/2l  »  ''^^ j-l/2^°' j-1/2   *    ^^^    » 

s     =     sgn(aH_^^/2)    5  (5.5c) 

Y^+1/2   =A^i/2g''   /  «5+i/2    ;  (5.5d) 

^j+1/2   =^   aj+1/2    '  ^^-^^^ 

here,  ct'^^w2  is  (5.4),  and  as  before,  Q(x)  is  (3.16)  and  a(x)  is 
(4.2b). 

lemma  5.1.  The  finite-difference  scheme  (5.5)  is  a  second  order 
accurate  approximation   to    (1.1). 

Proof :    For  each  k,    exactly  as   in   the    scalar  case,   we  have 


8j  +  gj+1   =  2  a(v^i/2)  «j+i/2   +  ^(^2) 


and 


Q(v^^/2   +  Y  j+1/2)  "j+1/2   =  Q(^  j+1/2)   "j+1/2   +  0(^    ^ 


hence 

2 
gj+  gj+i   -  Q(v  3+1/2  +  Tj+1/2)  «j+l/2   =    [2a(v^i/2)-   Q(v  j+1/2)]  aj+l/2  +  ^(^    > 

=   2(n    -  j)    (Vj+1/2)^  "j+1/2   • 

Using    (5.4a)   we    get 

m 
J    I      [g)+  g)+i   -  Q(v^+i/2  +  Yj+l/2)«j+l/2l^j+l/2 


k=l 
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k=l 
=    (n-1)   X2    [A(Vj+i/2)]^  ^jfl/2  ^  +  0(^2) 


Thus,    the   numerical    flux    (5.5b)    satisfies 


K    ,  2/ .2 


^    ^j+1/2    =  [X    f   +  h(n    -   2-)   ^    (A-u^)]  j^.i/2   +  O(h^), 

which   is    the    same   as    (4.3b).      The  rest    of    the      proof      is      identical      to 
that   of  Lemma    4. 1. 

We  turn  now  to  consider  the  constant  coefficient  case  A(u)  =  A.  = 
const.  in  (1.1).  In  this  case  the  eigenvalues  {a  }  and  the 
eigenvectors  [l^}  and  {r  }  are  also  constant,  and  the  pure  IVP 
decouples  into  ra  scalar  constant  coefficient  problems  for  the 
characteristic  variables,    i.e.. 


w*^   =  Jl^  V  (5.6a) 

w^  +  aK^  =  0   .  (5.6b) 


The  total  variation  in  x  of  each  of  the  characteristic  variables 
is  diminishing  in  time.  We  show  now  that  the  same  is  true  for  (5.5), 
i.e.  the  finite-difference  scheme  also  decouples  into  ra  scalar  ones 
for  the  characteristic  variables.  To  see  that,  we  multiply 
(5.5a)-(5.5b)    from  the   left   ty   i^  and  note   that 


^^(f j  +  f j+i)    =  i^   A(v.   +  vj+i)    =  ai(wi   +  w]^^)    ; 
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for  each  of   the   characteristic   variables    (dropping    the    superscript),    W2 
get 

w^-^^   +  n   X(f5:;/2    -  f~5!l/2)    =  wT   -   (1-n)   X(f~5^i/2-  f~5_i/2)  (5.7a) 

where 

f  j+1/2    =  ^    f  j+1/2    =  2   ^^^"j  ^  ""j+l^    "*"  X    ^^J   ^  ^J+1 

-   Q(v^l/2   +  Yj+1/2)   aj+1/2]}  (5.7b) 

Since 


■j+1/2    -^jfl/2 


(5.7c) 


it  follows   that    (5.7)  with  (5. 5c)-(5. 5e)    is    identical    in      the      constant 
coefficient      case      to      the      second      order  accurate    scalar  TVD  scheme   of 

Theorem  4.3.      Hence  it    is   TVD  with  TV(w)    =  J^       l«-i+l/2l      ^"^      convergent 

J 
under  the   CFL    restriction   (4.9).      We   conclude  that 

Theorem   5.2.        The      second      order     accurate      scheme      (5.5)      in  the 

constant  coefficient    case   is   TVD  and  convergent,   where 

TV(v)   =  I      I       |a^^i/2i    .  (^-^^ 

j   k=l 

subject    to   the   CFL    restriction 


ma 


k  ^  ^k       .    I    /      1      ^  (5.9) 


We  note  that  the  particular  form  of  averaging  in  (5.3)  does  not 
enter  into  tVe  considerations  of  Theorem  5.2.  However,  if  we  require 
the  schane  (5.5)  for  m  =  1  to  be  identical  to  the  scalar  scheme  of 
Section  4,    then  we  have   to  choose    (5.3)    so  that   "^  ^^.W2   i"   (5.5e)    is    the 
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same    as    the    mean   value    CFL   nunber    (3.10).      This   can  be  accomplished      by 
taking    the  eigenvalues   a^+W2   and  the   eigenvectors   '^^+1/2   and    r^+W2    ^^ 
(5.5)    to  be    those   of    A.(v .  ,v^^.i ) ,      where      A(u,v)      is     Roe's     mean     value 
Jacobtan.      This  matrix   satisfies    (see    [14]) 
(i)         f(v)   -   f(u)    =  A(u,v)(v-u) 
(ii)      A(u,u)    =  A(u) 
(iii)   A(u,v)   has  real   eigenvalues   and   a  complete    set    of   eigenvectors. 

The  existence  of  a  mean  value  Jacobian  that  satisfies  (iii) 
globally  is  related  to  the  symrae trizabili ty  of  (1.1)  and  can  be 
constructed  via  the  entropy  function  (see  [4]  and  [8]).  Roe  in  [14] 
constructs  a  mean  value  Jacobian  for  the  Euler  equations  of  gas 
dynamics  of  the  form  A(u,v)  =  A(V(u,v)),  where  V(u,v)  (5.3)  is  some 
particular  average. 

The  numerical  flux  (5.5b)  approximates  the  flux  across  x  =  0  in 
the  Riemann  problems  (5.2)  and  (5.4);  therefore,  the  scheme  (5.5)  can 
be  interpreted  as  a  second  order  accurate  Godunov-type  scheme  (see 
[3]). 
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6.    Applications   to   steady-state  calculations. 

Steady-state  solutions  (satisfying  f^  =  0)  for  the  initial  value 
problem  of  (1.1)  are  either  a  constant  state  (for  the  periodic  case)  or 
a  stationary  shock  (for  the  infinite  domain).  Hence  the  practical 
steady-state  problem  involves  a- forcing  term  (i.e.  f(^)x  ~  H(u,x))  and 
appropriate  boundary  conditions  to  maintain  a  richer  structure  of  a 
steady-state   solution. 

In  this  section  we  consider  the  calculation  of  a  steady-state 
solution  through  a  time -consistent  approach.  We  choose  initial  data 
u(x,0)    =  Uq(x)   such   that    u(x,t)    satisfies 

lim   u(x,t)    =  u„(x)    ,  (6.1) 

t-H» 

where  i^(x)  is  the  steady  state  solution,  and  advance  it  in  time  using 
a  scheme  of  the  form  (3.5a);  the  calculation  is  terminated  when 
llyn+l    _  vHii       falls      below      some      specified        tolerance.  Hence        this 

time-consistent  approach  involves  two  separate  limiting  processes:  (i) 
the  limit  n  -»■  "  with  fixed  h  and  t,  which  is  the  stationary  solution  of 
the  finite  difference  scheme,  lim  v"  =  Vj,  .  (11)  The  limit  h  -»•  0  of 
v^  ,  which  is  the  steady-state  solution  u„(x)  of  the  partial 
differential   equation. 

Indulging  in  linear  thinking,  we  may  consider  the  time-consistent 
approach  to  steady-state  to  be  computationally  equivalent  to  the  decay 
to  zero  of  the  solution  of  (1.1)  to  the  initial  data  u(x,0)  =  u^(x)  - 
Uq(x),  where  Ux,(x)  is  the  steady-state  solution  (6.1).  From  this  point 
of      view,      the      use      of      TVD     schemes      for      such      calculations      is   most 


-42- 
appropriate,    as  these      schemes      are      particularly     suitable      to      handle 
processes  of    decay. 

The  rate  of  decay  associated  with  the  pure  initial  value  problem 
(1.1)  is  0(t~^/2)  j^ri  the  infinite  domain,  and  0(t~^)  in  the  periodic 
case  (see.  [2]  and  [10]).  Hence  the  use  of  an  explicit  scheme  in  a 
time-consistent  approach  to  steady  state  may  turn  out  to  be 
prohibitively  expensive.  Therefore,  it  makes  good  computational  sense 
to  use  the  implicit  backward  Euler  scheme,  n  =  1  in  (5.5),  which  is 
unconditinally  TVD  (see  (5.9)).  This  enables  us  to  take  very  large 
time-steps  (that  are  guaranteed  not  to  increase  the  oscillation  of  the 
numerical  solution)  and  thus  obtain  fast  convergence  to  the  stationary 
solution  v^  of  the  finite  difference  scheme.  The  spatial  second  order 
accuracy  allows  v^  to  be  a  highly  resolved  yet  nonoscil  latory 
approximation    to    the    steady  state    solution   u^(x). 

We  turn  now  to  consider  the  second  order  accurate,  TVD,  modified 
backward  Euler  scheme    in    the    scalar  case: 

vfl   +UF5:;/2   -?5ll/2)    =  v]    ;  (6.2a) 


y-,_L.l 

^i+1/2      •'-^     ^^^      following     numerical    flux  calculated   at   the    (n+1)    time 


level 


^j+1/2    =  2    ^^J  ^  ^J+1    "^  2"  f^j  ^  "j+1    "   ^^''j+1/2   +  ^  j+1/2)  ^  j+1/2^]}  ;(6.2b) 


Q(x)    is   (3.16), 


J4-1/2V   =   Vj+i    -  vj    ; 
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^  j+1/2  =  ^  ^  j+i/2f/A  jfi/2v  ; 

^j+1/2    =  ^  jfl/28/^j+l/2V    ; 
gj    =  s    •    max  {O,    min[a(v  j^j/2)  1^  jf-1/2^1 '^•^'^^  j-1/2)  '^j-1/2^1}         (6.2c) 


and 


a(v) 


—  [Q(v)+v    ]    ,       for   time-dependent    problems 
,       for   steady-state   problems. 


(6. 2d) 


jQ(v) 


It  follows  from  (4. 9b)-(4. 9c)  that  the  scheme  (6.2)  can  be   expressed, 
as  in  (3.11),  by 

CT.-i,n+lN   -  ,,n+l  4.  /f+     x(n+l)A      ,,n+l  _  (r~  ^(n+l)  *       ,,ri+l  -  m^   ( f,    'ta^ 

(L-v    )j  -  Vj    +  (Cj^^/2)     ^j+1/2^      '■^j-1/2)      ^j-1/2  ^    -VjCb.Ja) 

where  by  (3.15b)  and  (4.9b) 


C±  =  -  C*  =i.  r±  (v  +  y)  -  Q(v  +  Y)] 


(6.3b) 


We  note  that  since  Q(x)  in  (3.16)  is  such  that  Q(x)  _>  |x|,  therefore  C~ 
in  (6.3b)  satisfy 


C?  <  0 


(6.4) 


If  we  take  Q(x)  =  |x|  (i.e.  e  =  0  in  (3.16))  then  C^  =  ±  (u  + 
Y  )"*" ,  which  corresponds  to  upstream  differencing  with  respect  to  the 
characteristic  field  (v  +y)'  (^n  the  right  hand  side  of  the  above 
equality  we  have  used  the  funcational  convention  b  =  max(b,0),  b  = 
min(b,0);  not  to  be  confused  with  the  notation  C-   on   the   left   hand 
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side.  )      Hence     Q(x)      with     e      >^     0     i  n      (6.2)      corresponds      to   upstream 
differencing    supplemented  by  a    viscosity   term  with   coefficient   e. 

The  backward  Euler  scheme  (6.3)  is  the  only  one  in  the  class  of 
schemes  (3.5)  to  have  the  identity  operation  as  its  right  hand  side. 
Since  it  is  the  requirement  for  the  right  hand  side  to  be  a  TVD 
operator  that  limits  the  CFL  nunber,  it  is  also  the  only  scheme  in  the 
class   (3.5)    that  is    unconditionally  TVD. 

The  scheme  (6.2)  is  implicit,  i.e.,  the  value  of  v""*"^  is  obtained 
as  the  solution  of  a  system  of  nonlinear  algebraic  equations.  This 
system  consists  of  the  equations  (6.2)  for  all  j  in  the  computational 
domain,  supplemented  by  appropriate  boundary  conditions.  (To  simplify 
our  presentation  we  assume  periodic  boundary  conditions  throughout  this 
section.)  The  iteration  needed  to  solve  this  system  of  nonlinear 
equations,  at  each  time-step,  say,  by  Newton's  method  may  completely 
wipe   out    the   savings    gained  by  the  ability    to   take    larger   time-steps. 

To  overcome  this  obstacle,  we  consider  a  linearized  version  of 
(6.3)  in  which  the  coefficients  (C-)^""'"!)  are  replaced  by  (C*)^"\ 
i.e.  , 

(L.vn+l)j  ,    vf  1    +    (Ct^i/2)^"^   Vl/2-"^'    "    ^^j-l/z)^"^   ^j-1/2   v"^^    =  v'j    (6.5a) 

or 

,,x      _    (r-  x(n),,n+l    4.   n    -   r~  p+  x(n),,n+l    .    (-p+      ,    \(n)..n+l    =  .,n 

(6.5b) 

(6.5)  is  a  tridiagonal  system  of  linear  equations;  furthermore,  the 
nonpositivi ty  of  the  coefficients  C-  (6.4)  implies  that  (6.5b)  is 
diagonal  ly  dominant .      Therefore    it    follows   immediately   that 
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LemiTH  6.  1.  L,  the  finite  difference  operator  on  the  left  nand 
side   of    (6.5)    is  unconditionally  invertible. 

The  nonpositi vi ty  of  C"  in  (6.4)  also  implies  that  the  operator  L 
in  (6.5)  is  TVI  (see  part  (b)  of  LertiTiH  3.2).  Since  the  right  hand  side 
of  (6.5)  is  the  identity  operator,  it  follows  immediately  from  Lemma 
3.1    that 

Lemna  6.2.  The  linearized  backward  Euler  implicit  scheme  (6.5)  is 
unconditionally  TVD. 

We  note  that  the  linearized  scheme  (6.5)  is  not  in  conservation 
form  and  therefore  should  not  be  used  to  approximate  time-dependent 
solutions.  The  following  theorem  shows  that  it  is,  however,  a  suitable 
scheme    for  the    calculation   of  steady-state  solutions. 

Theorem  6.3.  If  v^^  is  the  limit  n  ^-  °° ,  with  h  and  t  fixed,  of 
(6.5),    i.e.      lira  v"  =   v^    ,    then   v*  satisfies 

fj+l/2(v*)    -   f j-i/2(v*)    =  0  (6.6) 

where    ^i+i/2(^*^   is    the   numerical    flux    (6.2b)    computed   at   v*    . 
Proof :    Rewrite    (6.5a)    in    the    form 

^i  +  (ct^l/2)^"^   ^j+l/2d  -    (C]-i/2)^"^   A  ._i/2   d 

=   (C+^i/2)^"^   Aj+1/2-"-    (CT_i/2)^"^   Aj-1/2   v"  (6.7a) 


where 


d.   =   v^+l    -   v^  (6.7b) 
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and  observe    that    by    (4.9b)    the   right  hand    side  of    (6.7a)      is     equal      to 
-  ■^(f  j+1/2"   fj-l/2)>    whei^e    fj+i/2    is   (6.2b),    i.e. 

dj+    (ct^i/2)^"^   A^l/2  d  -    (C]-i/2)^"^   A._i/2   d 

=   -A(?5+i/2   -i'^-i/2).  (6.7c) 

Next        we      note      that      v*^      is      a      finite-dimensional      vector,      and 

therefore,    the    convergence    v'^  +    v*   is    componentwise.      Hence    lim   d^    =      0 

n-H» 
for  all    j   and    (6.6)    follows   from    (6.7c). 

We      turn     now     to     describe      the      extension  of    (6.5)    to  systems   of 

coraservation   laws.      As   in  Section    5,    our   only  design  guideline      is      the 

requirement      that      the      scheme      will  be  TVD  in  the   constant    coefficient 

system   case.      To    ensure   that    the   steady  state      solution      is      consistent 

with      the     conservation      form      in     the      sense   of    (6.6),    we  use    the    form 

(6.7c)    rather   than   (6.5).      Using    the   extension   technique  of    Section      5, 

we    get 

"^j  "^  ^j+l/2  A  j+i/2  d  -    l^j-1/2  Aj-l/2   d   =    -AUj+i/2-f'j-l/2)  (6.8a) 

where    the  numerical    flux   f'^  is    (5. 5b)-(5. 5c)  with  0(x)    given  by   (6. 2d); 
Kr   are   rmm  matrices   given  by 

+  +  —  1 

^"j+1/2   =  Rj+l/2  'l^l+l/Z   Rj+1/2  (6.8b) 

where   R-;+W2   =   R(vi+i/2)    (5.1a)    and  \l)~    are    the   diagonal  matrices 
('^j+l/2>i,k  =  y[±(^j+1/2   +vVl/2)    -  Q(v^+1/2+yVi/2)]    «i,k   .    (6.8c) 
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The  scheme    (6.8)   can  also  be   written    In   the   operator   form 

(I  +  K+^i/2  A^i/2   -  KT_i/2  ^^  j-1/2)  (^""'^    -  v")    =   -  X(F]+i/2   "  f"j-l/2)    (6.9) 

where   I-w^,  Is   the   identity  matrix. 

Next  we  consider  the  constant  coefficient  case  and  note  that  (6.9) 
differs  from  (5.5)  with  n  =  !•  IVie  reason  is  that  (5.5)  is  nonlinear 
even  in  the  constant  coefficient  case,  and  (6.9)  is  its  linearized 
version.  However,  since  (6.9)  decouples  into  m  scalar  schemes  of  the 
form  (6.5)  for  the  characteristic  variables,  we  conclude  by  Lemma  6.2 
and  Theorem  2.1   that 

Theorem  6.4.  The        second        order       accurate        implicit      scheme 

(6.8)-(6.9)  in  the  constant  coefficient  case  is  unconditinally  TVD 
(5.8)   and   convergent. 

Numerical  experiments  in  [17]  demonstrate  the  superiority  of  (6.9) 
over  conventional  methods. 

Summary 

A  conventional  scheme  for  the  solution  of  nonlinear  hyperbolic 
conservation  laws  is  one  that  is  linear  and  L2-stable  when  considered 
in  the  constant  coefficient  case.  A  typical  complaint  about  the 
application  of  conventional  schemes  to  the  solution  of  nonlinear 
problems    is  that   they  are  either  not    robust   and/or  not  accurate  enough. 

In  this  paper  we  propose  a  more  appropriate  setting  for  the  design 
and  analysis  of  schemes  for  nonlinear  problems,  in  which  one  can  obtain 
both   robustness   and   accuracy  at  the   same   time. 

We  use   the      convergence      Theorem     2.1      to     rigorously      settle      the 


question  of  "robustness".  Under  the  guidelines  of  this  theorem,  we 
consider  the  class  of  schemes  that  are  total-variation  stable  and 
consistent   with   the   conservation    law   and    its    entropy    inequality. 

Vfe  present  a  technique  to  convert  first  order  accurate  schemes  in 
this  class  into  second  order  ones.  We  remark  that  2  seems  to  be  the 
maximal  order  of  accuracy  of  total-variation  diminishing  schemes,  since 
monotonicity  requirements  seem  to  necessitate  the  occasional 
deterioration  of   the    truncation  error    to   OCh"'). 

The  theory  of  this  class  of  schemes  is  still  in  its  infancy.  Much 
more  work,  is  needed  to  formulate  necessary  and  sufficient  conditions 
for  uniform  boundedness  of  the  total-variation,  as  well  as  for 
consistency  with   the  entropy    inequality. 
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Appendlx 
by   Peter  D.   Lax 

We  give  here   a    potentially  useful  general   criterion    for      a      linear 

difference      operator      to      diminish      total   variation.      We  start  with   the 

general   observation   that    over  Z,    the    i      and  i      norms    are   dual  to      each 

other  with   respect    to   the   Jl^  duality.      It    follows   then   that    if  S    is  any 

operator  mapping    functions    on  Z   into    functions    on  Z,    and    if   S  is      the 

transpose      of      S      with      respect      to      the      i^      scalar      product,  then   the 
following    relation   holds  between   the   norms   of    S   and   S    : 


isl,,  -  ls*l^. 


We  take  now  S  to  be  a  difference  operator; 


(1)  S  =  ^  a^TJ  , 


T    translation  on  Z ,    a-   functions    on  Z.    Then 


(D*  S*   =  I   T-Ja.   =  I   a^jj^  T-J    , 

where 

(2)  aU)(k)    =  aj(k  +   j) 

Now    tie   l°°  norm  of   a  difference   operator  is    easily   determined; 


(3)  1S*|   „    =Max  );    |aU\k) 


'  k      j 


so    |S*|         <   1   iff 
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(4)  I     |aj(k  +  j)  I    <   1      for  all     k. 

J 

This   is    in   particular   the    case    if 

(4)'  3j   >    0     and     I    aj(k+j)    <   1      for  all  k   . 

j 

Thus,    (4)    is  the   necessary  and    sufficient    condition    for  S    given  by      (1) 
to   diminish    the  i^   norm. 

Next      we      show  how  to  use    this    to  deduce  a    condition    that    a    linear 
difference  operator 


(5)  R  =  I    bjTJ 


diminish    total   variation.      By   definition,    for  any  function   v  defined  on 
Z, 


(6)  TV   V  =       I (I   -    T)v      , 

I'- 


We     suppose      now     that   R  preserves   the   constant    function.      This   is 
the   case   iff 

(7)  I   bj(k)    =   1     for  all  k. 

j 

Then  there    is    an  operator  S   such   that 

(8)  (I-T)R  =  S(I    -  T) 

To    see   this    we    express   both  sides   as    difference   operators: 
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(I-T)R  =   y    b,TJ    -  y    T   b-TJ 

=1    bjTJ    -I    b(-l)TJ+l       =     I    (bj    -b(:}))T3. 
Similarly 

S(I-T)      =     I   a^TJ   -  a^TJ+l    =1    (a^   -  aj_i)TJ    . 
Therefore    (8)    holds  iff 


( 


I.e. 


(9) 


aj(k)    -  aj_i(k)      =     bj(k)    -  bj_i(k-l) 


From  these  equations  we  can  determine  a. 
(10)    aj(k)   =  I   b^Ck)  -  I      b.(k-l). 


We  claim  that  if  the  operator  R  is  explicit,  i.e.  if  b-  =  0  for  |j|  > 
J,  then  so  is  the  operator  S.  For  formula  (10)  shows  that  then  a-  =  0 
for  j  <  -J  ;  combining  (10)  with  (7)  we  see  that  a.  =  0  also  for  j  >  J. 

The  same  argument  shows  that  if  R  is  implicit,  i.e.  b.  ^  0  as  j  > 
±<*> ,    the  same  is  true  for  the  operator  R. 

We  have  shown: 

The  operator  R  given  by  (5),  satisfying  condition  (7),  is  total 
variation  diminishing  iff  the  operator  S,  whose  coefficients  are  given 
by   (10),    satisfies    condition   (4). 
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